When an SDEProblem is given a sparse noise_rate_prototype, _sde_init built the Wiener increment prototype as false .* noise_rate_prototype[1, :], which is a SparseVector. That is the wrong shape for the object: the sparsity pattern of g records which states a channel drives, not which channels are idle, so every column still draws an increment and dW fills in completely on the first step.

Most solvers cope with this, but two do not.

W2Ito1 fails outright. Its cache keeps _dW, _dZ and chi1 under a single randType field, and builds _dZ = zeros(eltype(ΔW), 2) unconditionally. With a sparse dW the three arguments have two different types and there is no applicable constructor:

using StochasticDiffEq, SparseArrays

f!(du, u, p, t) = (du .= -u)
function g!(du, u, p, t)
    du[1, 1] = 0.1 * u[1]
    du[3, 2] = 0.1 * u[3]
end

np = spzeros(4, 2); np[1, 1] = 1.0; np[3, 2] = 1.0
prob = SDEProblem(f!, g!, ones(4), (0.0, 1.0), noise_rate_prototype = np)

solve(prob, W2Ito1(), dt = 0.01, adaptive = false)

ERROR: MethodError: no method matching StochasticDiffEqWeak.W2Ito1Cache(
  ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64},
  ::SparseVector{Float64, Int64}, ::Vector{Float64},
  ::SparseVector{Float64, Int64}, ...)

RKMilGeneral fails silently, which is worse. It derives its Lévy area dZ from the increment prototype through similar, so a sparse prototype gives a sparse dZ of length 43 for p = 10 and two noise channels. randn! has a specialized method for Array{Float64} and a generic AbstractArray loop for everything else, and past a small length the two consume the RNG in a different order. Same problem, same seed, sparse versus dense noise_rate_prototype:

sparse prototype: u[end][1] = 0.34044131290011104
dense  prototype: u[end][1] = 0.34940768585718085

The fix folds the sparse noise_rate_prototype case into the existing issparse(u) branch, which already builds a dense increment through adapt(SciMLBase.parameterless_type(u), zeros(...)). Reusing that branch rather than a bare zeros keeps the array type of u, so a GPU u still gets a device-resident increment.

Keeping the increment dense is also cheaper, because the increment type is what the saved noise history stores. 400 states over 200 channels, EM at dt = 1e-3 out to t = 10:

master: 136.71 MiB, 58 ms
branch:  33.25 MiB, 25 ms

Testing

New file lib/StochasticDiffEq/test/sparse_noise_tests.jl, registered in runtests.jl. Run against unmodified master it gives 17 passed, 3 failed, 1 errored; on this branch 21 passed. The four that move are the ones above: integrator.W.dW is a Vector{Float64} rather than a SparseVector, W2Ito1 reaches ReturnCode.Success instead of throwing, and RKMilGeneral(p = 10) produces the same trajectory from a sparse prototype as from Matrix(noise_prototype) with the same seed. The remaining testsets cover EM and EulerHeun in place and out of place and pass either way.

Also ran the nearby suites on the branch with no change: nondiag_noise_eulerheun_test 2/2, sparsediff_tests 5/5, nondiagonal_tests 5 passed 2 broken, noise_type_test 2024/2024, scalar_noise 1/1.

This addresses part of #3178. That issue also asks for sparse noise processes that only generate the channels a problem actually uses, which this does not do, so I have not marked it closed.

Two adjacent instances of the same false .* noise_rate_prototype[1, :] pattern are left alone, both in DelayDiffEq: ext/DelayDiffEqStochasticDiffEqCoreExt.jl:48 and src/utils.jl:435 (dW_dummy). They have the identical defect but sit in a different sublibrary, so they are out of scope here.

AI Disclosure

Claude was used to help investigate the failure and draft this change.